program test

    use iasing_osse
    use netcdf

    implicit none
    integer :: i, jl, jo
    ! netcdf vars
    character(len=128)  :: ncfile, cnopro
    integer(JPI4)       :: ncid, nopro, nlevels, nchl
    integer             :: chia, chng
    type(r4a_ids)       :: db_ids
    type(r51_dat)       :: db_dat
    real(JPR4)          :: tia(NIASI), tngA(NIANG), tngB(NIANG), emis(NIANG)

    integer :: nargs, nobs
    integer, allocatable :: obslist(:)
    character(len=128) :: cargs(10), list
    logical :: skip_next = .false.
    logical :: list_file = .false.
    logical :: brt_units = .false.

    integer :: jout ! output unit

    include "checknc.interface"
    include "r51_initializer.interface"
    include "r51_reader.interface"
    include "rad2brt_f.interface"
    include "brt2rad_f.interface"

    data jout /6/

    nargs = command_argument_count()
    if (nargs.eq.0) then
        print*, "No arguments supplied. Aborting execution..."
        stop -1
    else
        ! Read arguments
        do i = 1,  nargs
            call get_command_argument(i, cargs(i))
            !print*, i, trim(cargs(i))
        enddo
        ! Process arguments
        do i = 1, nargs
            if (skip_next .eqv. .true.) then
                skip_next = .false.
                cycle
            endif
            if (cargs(i).eq."-f") then
                skip_next = .true.
                read(cargs(i+1), '(A)') ncfile
            else if (cargs(i).eq."-o") then
                skip_next = .true.
                read(cargs(i+1), *) nopro
            else if (cargs(i).eq."-l") then
                read(cargs(i+1), '(A)') list
                skip_next = .true.
                list_file = .true.
            else if (cargs(i).eq."-u") then
                skip_next = .true.
                read(cargs(i+1), *) jout
            else
                print*, "Error in the supplied arguments. Please check the doc..."
                print*, "Aborting execution"
                stop -2
            endif
        end do
    endif

    ! Read list of observations, just in the case it has been provided
    if (list_file .eqv. .true.) then
        open(1, file=list, action='read', form='formatted', status='old')
        read(1,'(A)') ncfile
        read(1,*) nobs
        if (nobs.gt.0) then
            allocate(obslist(nobs))
            read(1,*) obslist
        else
            call checknc(nf90_inquire_dimension(ncid, db_ids%atm_did, len=nobs))
            print*, nobs
            allocate(obslist(nobs))
            obslist = (/(i,i=1, nobs)/)
        endif
        close(1)
    else
        nobs = 1
        allocate(obslist(nobs))
        if (nopro/=0) then
            obslist = nopro
        else
            obslist = 1
        endif
    endif

    ! Open netcdf file
    call checknc(nf90_open(ncfile, NF90_NOWRITE, ncid))
    db_ids = r51_initializer(ncid)
    call checknc(nf90_inquire_dimension(ncid, db_ids%lev_did, len=nlevels))
    call checknc(nf90_inquire_dimension(ncid, db_ids%ocl_did, len=nchl))

    write(*,*)
    write(*,'(2A)') "ncfile: ", trim(ncfile)
    write(*,*)

    do jo = 1, nobs
        ! Read record from database
        nopro = obslist(jo)
        db_dat = r51_reader(ncid, db_ids, nopro)
        
        ! Obtain brightness temperatures and emissivity
        emis = real(db_dat%opa%ert,4)*.001

        ! Change Get Brightness temperature and change units when necessary
        ! IASI
        if (any(db_dat%rbt%ria<0. .or. db_dat%rbt%ria>1000.)) then
            db_dat%rbt%ria = -999.9
            tia = -999.9
        else if (db_dat%rbt%ria(1)>200.) then
            tia  = db_dat%rbt%ria
            db_dat%rbt%ria = brt2rad_f(NIASI, db_dat%dim%wia, tia)
        else
            tia = rad2brt_f(NIASI, db_dat%dim%wia, db_dat%rbt%ria)
        endif
        ! IASI-NG A
        if (any(db_dat%rbt%rng<0. .or. db_dat%rbt%rng>1000.)) then
            db_dat%rbt%rng = -999.9
            tngA = -999.9
        else if (db_dat%rbt%rng(1)>200.) then
            tngA  = db_dat%rbt%rng
            db_dat%rbt%rng = brt2rad_f(NIANG, db_dat%dim%wng, tngA)
        else
            tngA = rad2brt_f(NIANG, db_dat%dim%wng, db_dat%rbt%rng)
        endif
        ! IASI-NG B
        if (any(db_dat%rbt%ngb<0. .or. db_dat%rbt%ngb>1000.)) then
            db_dat%rbt%ngb = -999.9
            tngB = -999.9
        else if (db_dat%rbt%ngb(1)>200.) then
            tngB  = db_dat%rbt%ngb
            db_dat%rbt%ngb = brt2rad_f(NIASI, db_dat%dim%wng, tngB)
        else 
            tngB = rad2brt_f(NIANG, db_dat%dim%wng, db_dat%rbt%ngb)
        endif

        ! Write results
        write(jout,'("#",2(A,I8))') "Profile number: ", db_dat%dim%nopro, &
            ", corresponding index-list: ", jo
        write(jout,'("# ",A)') " COORDINATES "
        write(jout,'("# ",2(A,I10))') "Date: ", db_dat%coo%dat, ", time: ", &
            db_dat%coo%tim
        write(jout,'("# ",3(A,F9.3))') "Longitude: ", db_dat%coo%lon, &
            ", latitude: ", db_dat%coo%lat, ", elevation: ", db_dat%coo%ele/G0
        write(jout,'("# ",A)') " OBSERVATION PARAMETERS "
        write(jout,'("# ",2(A,F9.3))') "Surface pressure: ", db_dat%opa%psu, &
           ", Surface temperature: ", db_dat%opa%tsu
        write(jout,'("# ",2(A,F9.3))') "Satellite, Zenith angle: ", &
            db_dat%opa%iza, ", Azimuth angle: ", db_dat%opa%iaa
        write(jout,'("# ",2(A,F9.3))') "Solar,     Zenith angle: ", &
            db_dat%opa%sza, ", Azimuth angle: ", db_dat%opa%saa
        write(jout,'("# ",A,F8.4,A,I6)') "Land-sea mask: ", &
            db_dat%opa%lsm, ", Cloud cover: ", db_dat%opa%clc
        write(jout,'("# ",A)') " VERTICAL PROFILES "
        write(jout,'("# ",A7,A10,6A12)') "P", "T", "q", "co2", "o3", "co", &
            "ch4", "so2"
        do jl = 1, nlevels
            write(jout,'(F9.4,F10.3,6ES12.4)') db_dat%dim%lev(jl), &
                db_dat%vpr%tem(jl), db_dat%vpr%q(jl), db_dat%vpr%co2(jl), &
                db_dat%vpr%o3(jl), db_dat%vpr%co(jl), db_dat%vpr%ch4(jl), &
                db_dat%vpr%so2(jl)
        enddo
        write(jout,'("# ",A)') " CLOUD INFORMATION "
        write(jout,'("# ",A7,5A12)') "P", "cc", "iwc", "lwc", "rwc", "swc"
        do jl = 1, nlevels
            write(jout,'(F9.4,5ES12.4)') db_dat%dim%lev(jl), &
                db_dat%cld%cc(jl),  db_dat%cld%iwc(jl), db_dat%cld%lwc(jl), &
                db_dat%cld%rwc(jl), db_dat%cld%swc(jl)
        enddo
        write(jout,'("# ",A)') " IASI OBSERVATIONS "
        write(jout,'("# ",A7,6A10)') "Channel", "obs", "fgd", "bia", "iasi", &
            "iang_A", "iang_B"
        do jl = 1, nchl
            chia = db_dat%dim%ocl(jl)
            chng = chia*2-1
            write(jout,'(I9,3F10.3,3F10.3)') db_dat%dim%ocl(jl), &
                db_dat%obs%obt(jl), db_dat%obs%ofd(jl), db_dat%obs%obi(jl), &
                tia(chia), tngA(chng), tngB(chng)
        enddo
        write(jout,'("# ",A)') " SIMULATIONS "
        write(jout,'("# ",A6,A8,3A12,3A10)') "wnu", "emis", "ria", "rngA", "rngB", &
            "tia", "tngA", "tngB"
        do jl = 1, NIANG
            if (mod(jl,2).eq.1) then
                chia = jl/2+1
                write(jout,'(2F8.3,3ES12.4,3F10.3)') db_dat%dim%wng(jl), &
                    emis(jl), db_dat%rbt%ria(chia), db_dat%rbt%rng(jl), &
                    db_dat%rbt%ngb(jl), tia(chia), tngA(jl), tngB(jl)
            else
                write(jout,'(2F8.3,3ES12.4,3F10.3)') db_dat%dim%wng(jl), &
                    emis(jl), -999.9, db_dat%rbt%rng(jl), db_dat%rbt%ngb(jl), &
                    -999.9, tngA(jl), tngB(jl)
            endif
        enddo



        ! Radiances
        !print*
        !print '(A)', "Instrument spectra at IASI wavenumbers"
        !print '(A,10F10.3)', "ria(1:10)", db_dat%rbt%ria(1:10)
        !print '(A,10F10.3)', "rng(1:10)", db_dat%rbt%rng(1:20:2)
        !print '(A,10F10.3)', "ngb(1:10)", db_dat%rbt%ngb(1:20:2)

        write(jout,*)
        write(jout,*)
    enddo

    call checknc(nf90_close(ncid))

end program
